import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
import sys
import os

# 设置中文字体
import matplotlib
matplotlib.rcParams['font.family'] = ['SimHei']
matplotlib.rcParams['axes.unicode_minus'] = False

# 创建输出目录
os.makedirs('C', exist_ok=True)

# 设置输出重定向
LOG_FILE = 'C/C.txt'
log_file = open(LOG_FILE, 'w', encoding='utf-8')
original_stdout = sys.stdout
sys.stdout = log_file

print("="*100)
print("古代玻璃制品成分分析与鉴别 - 问题3：未知玻璃样本分类分析")
print("="*100)

# 1. 数据加载
print("\n步骤1: 数据加载")
print("-"*100)

# 加载统计数据
stats_sheets = pd.read_excel("handled/四类文物统计.xlsx", sheet_name=None)
unknown_samples = pd.read_excel("附件.xlsx", sheet_name="表单3")

print(f"未知样本数量: {len(unknown_samples)}")

# 定义化学成分列表
composition_cols = [
    "二氧化硅(SiO2)", "氧化钠(Na2O)", "氧化钾(K2O)", "氧化钙(CaO)",
    "氧化镁(MgO)", "氧化铝(Al2O3)", "氧化铁(Fe2O3)",
    "氧化铜(CuO)", "氧化铅(PbO)", "氧化钡(BaO)", "五氧化二磷(P2O5)",
    "氧化锶(SrO)", "氧化锡(SnO2)", "二氧化硫(SO2)"
]

# 2. 主分类规则提取
print("\n步骤2: 主分类规则提取")
print("-"*100)

def classify_glass(row):
    """基于化学成分的玻璃类型分类"""
    if row['氧化铅(PbO)'] <= 5.46:  # 从B.py中得到的阈值
        return '高钾'
    else:
        return '铅钡'

# 3. 数据预处理
print("\n步骤3: 数据预处理")
print("-"*100)

# 处理未知样本数据
unknown_samples[composition_cols] = unknown_samples[composition_cols].fillna(0)

# 计算成分总和
unknown_samples['成分总和'] = unknown_samples[composition_cols].sum(axis=1)
valid_samples = unknown_samples[
    (unknown_samples['成分总和'] >= 85) & 
    (unknown_samples['成分总和'] <= 105)
].copy()

print(f"有效样本数: {len(valid_samples)}/{len(unknown_samples)}")

# 4. 样本分类
print("\n步骤4: 样本分类")
print("-"*100)

valid_samples['预测类型'] = valid_samples.apply(classify_glass, axis=1)
print("\n分类结果统计:")
print(valid_samples['预测类型'].value_counts())

# 5. 分类结果验证
print("\n步骤5: 分类结果验证")
print("-"*100)

def validate_classification(sample, class_stats_df):
    """验证样本分类结果的可靠性"""
    validation_scores = {}
    
    for comp in composition_cols:
        sample_value = sample[comp]
        # 获取当前成分的统计值
        stats_row = class_stats_df[class_stats_df['成分'] == comp]
        
        if not stats_row.empty:
            mean_value = stats_row['均值'].iloc[0]
            var_value = stats_row['方差'].iloc[0]
            
            # 计算z-score
            if var_value > 0:  # 避免除以零
                z_score = (sample_value - mean_value) / var_value**0.5
                validation_scores[comp] = abs(z_score)
    
    return np.mean(list(validation_scores.values())) if validation_scores else np.nan




# # 在主代码中添加调试信息
# print("\n数据结构调试:")
# for sheet_name in stats_sheets:
#     print(f"\n{sheet_name} 表的数据结构:")
#     print(stats_sheets[sheet_name].head())
#     print("\n列名:", stats_sheets[sheet_name].columns.tolist())


# 对每个样本进行验证
validation_results = []
for idx, sample in valid_samples.iterrows():
    pred_type = sample['预测类型']
    weathering = sample['表面风化']
    sheet_name = f"{pred_type}{'风化' if weathering == '风化' else '无风化'}"


    
    if sheet_name in stats_sheets:

        validation_score = validate_classification(
            sample, 
            stats_sheets[sheet_name]
        )
    else:
        validation_score = np.nan
        
    validation_results.append({
        '样本编号': sample['文物编号'],
        '预测类型': pred_type,
        '风化状态': weathering,
        '可靠性得分': validation_score
    })

validation_df = pd.DataFrame(validation_results)
print("\n分类可靠性统计:")
print(validation_df.groupby(['预测类型', '风化状态'])['可靠性得分'].describe())

# 6. 敏感性分析
print("\n步骤6: 敏感性分析")
print("-"*100)

def sensitivity_analysis(sample, threshold_range):
    """对分类阈值进行敏感性分析"""
    original_class = classify_glass(sample)
    results = []
    
    for threshold in threshold_range:
        temp_sample = sample.copy()
        temp_sample['氧化铅(PbO)'] = threshold
        new_class = classify_glass(temp_sample)
        
        results.append({
            'threshold': threshold,  # 注意这里使用英文列名
            'changed': original_class != new_class,
            'original_class': original_class,
            'new_class': new_class
        })
    
    return pd.DataFrame(results)

# 对关键样本进行敏感性分析
threshold_range = np.linspace(4.5, 6.5, 21)
sensitivity_results = []

# 分析每个样本
for idx, sample in valid_samples.iterrows():
    if 4.5 <= sample['氧化铅(PbO)'] <= 6.5:  # 临界样本
        result_df = sensitivity_analysis(sample, threshold_range)
        result_df['sample_id'] = sample['文物编号']
        sensitivity_results.append(result_df)

# 合并所有结果
if sensitivity_results:
    sensitivity_df = pd.concat(sensitivity_results, ignore_index=True)
    print("\n分类阈值敏感性分析:")
    print("阈值每变化0.1时的分类变化率:")
    print(sensitivity_df.groupby('threshold')['changed'].mean())
else:
    print("\n没有找到在阈值范围内的临界样本")

# 7. 结果可视化
print("\n步骤7: 结果可视化")
print("-"*100)

# 7.1 分类结果散点图
plt.figure(figsize=(12, 8))
sns.scatterplot(
    data=valid_samples,
    x='氧化铅(PbO)',
    y='氧化钾(K2O)',
    hue='预测类型',
    style='表面风化',
    s=100
)
plt.title('未知样本分类结果散点图')
plt.axvline(x=5.46, color='r', linestyle='--', label='分类阈值')
plt.legend(title='')
plt.grid(True, linestyle='--', alpha=0.7)
plt.savefig('C/classification_scatter.png')
plt.close()

# 8. 结果输出
print("\n步骤8: 结果输出")
print("-"*100)

# 准备输出结果
output_results = valid_samples[[
    '文物编号', '预测类型', '表面风化', '成分总和'
] + composition_cols].copy()

# 添加可靠性得分
output_results = output_results.merge(
    validation_df[['样本编号', '可靠性得分']],
    left_on='文物编号',
    right_on='样本编号',
    how='left'
)

# 保存结果
output_results.to_excel('C/classification_results.xlsx', index=False)
print("\n分类结果已保存至 C/classification_results.xlsx")

# 关闭输出重定向
sys.stdout = original_stdout
log_file.close()

print(f"分析完成！详细结果已保存至 {LOG_FILE}")
print(f"可视化结果保存至 C/ 目录")